Limits of validity for a semiclassical mean-field two-fluid model for Bose-Einstein 

condensation thermodynamics 
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We reinvestigate the Bose-Einstein condensation (BEC) thermodynamics of a weakly interacting 
dilute Bose gas under the action of a trap using a semiclassical two-fluid mean-field model in order 
to find the domain of applicability of the model. Such a model is expected to break down once the 
condition of diluteness and weak interaction is violated. We find that this breakdown happens for 
values of coupling and density near the present experimental scenario of BEC. With the increase of 
the interaction coupling and density the model may lead to unphysical results for thermodynamic 
observables. 
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I. INTRODUCTION 

The recent experimental observation of Bose-Einstein 
condensation (BEC) in a weakly interacting dilute gas of 
87 Rb 0, 23 Na §, 7 Li @, and X H f| employing magnetic 
traps at ultra-low temperatures calls for a theoretical in- 
vestigation on various aspects of the condensate. The 
condensate consists of few thousands to few millions of 
atoms confined by the trap potential. As the temperature 
is lowered below the critical temperature To of BEC, the 
condensate starts to form and finally at K all the avail- 
able atoms will be condensed. In the absence of a micro- 
scopic equation, the condensate is usually described by 
the mean-field Gross-Pitaevskii (CP) equation j5j. One 
of the primary interest on the process is to study the 
thermodynamic observables of the system, such as, the 
condensate fraction, internal energy, and specific heat. 

There have been several comprehensive studies on tem- 
perature dependencies of thermodynamic observables of 
the condensate using semiclassical mean-field two-fluid 
models The physical ingredients of these models 

are quite similar. One such model using the GP wave 
function provides satisfactory description of the temper- 
ature dependencies of the thermodynamic observables in 
two ||, and three |||l2| dimensions. These studies em- 
ployed a iterative solution of the system of equations in- 
volved. For the condensation of a system composed of 
40000 trapped 87 Rb atoms with repulsive interatomic in- 
teraction the iterative scheme converged rapidly and pro- 
vided a satisfactory account of the condensate fraction, 
internal energy, and specific heat ||. Similar conclu- 
sion was also reached in the study of condensation of 7 Li 
atoms JT^] . In case of 7 Li the attractive interatomic inter- 
action is responsible for collapse if the number of atoms 
is larger than approximately 1400 |D| . 

Here we reinvestigate critically the BEC of a weakly 
interacting dilute gas in two and three dimensions using 
the two-fluid mean-field model mentioned above in or- 
der to define the domain of its applicability. We have 
included the two-dimensional BEC in this study because 



of considerable recent interest in this topic |^[l4-17|. We 
employ the usual iterative solution of the nonlinear two- 
fluid mean-field model and study the convergence of the 
iterative scheme. Although, the convergence is rapid for a 
weakly interacting dilute system, with an increase of the 
strength of interaction and/or density, the model breaks 
down and leads to physically unacceptable results for the 
thermodynamic observables. Specifically, below the crit- 
ical temperature, the model may yield negative specific 
heat. 

However, it is well-known that the mean-field descrip- 
tion of the condensate via the GP equation as used in the 
semiclassical models above should hold under the condi- 
tion of diluteness of a weakly interacting Bose gas and 
is expected to break down once the conditions of dilute- 
ness and weak interaction are violated Q. The breakdown 
should happen for a large number of condensed particles 
reflecting a large density as well as for a large modulus 
of the scattering length denoting a strong interatomic in- 
teraction. These two conditions correspond to a large 
nonlinearity of the system and may lead to a breakdown 
of the mean-field two-fluid thermodynamic model p|,pT[ . 
In addition, the finite size of the system may necessitate 
corrections in the thermodynamic model as we are away 
from the real thermodynamic limit N — > oo, V — > oo, 
where N is the number of particles and V the volume of 
the system f\ Nevertheless exact numerical conditions for 
the breakdown of the mean-field thermodynamical model 
have never been investigated. In this work by performing 
numerical calculations we identify such conditions. We 
find that the semiclassical model may break down un- 
der possible present experimental conditions of BEC of a 
trapped Bose gas. 

We present the semiclassical model in Sec. 2, numerical 



*See, for example, page 474 of Ref. jLl| 

1 These limitations of the thermodynamic model are dis- 
cussed in Sec. 5 of Ref. 
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results in Sec. 3, and conclusions in Sec. 4. 



II. MEAN-FIELD MODEL 

We consider a system of iV bosons with repulsive in- 
teraction at temperature T under the influence of a trap 
potential. The condensate is described by the following 
GP wave function in the Thomas-Fermi approximation 
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) -90i-V ext -2gni{r)) i 
(2.1) 



where 9(x) is the step function, 6(x) = if x < and 1 
otherwise. Here V ex ^(r) = mw 2 r 2 /2 is the trap poten- 
tial, g the strength of the repulsive interaction between 
the atoms, m the mass of a single bosonic atom, lu the 
angular frequency, p the chemical potential, and n\{r) 
represents the distribution function of the noncondensed 
bosons. As we are interested in studying the limits of 
validity of the semiclassical mean-field model and not in 
simulating a particular experimental situation, we con- 
sider a spherically symmetric trap both in two and three 
dimensions. The noncondensed particles are treated as 
non-interacting bosons in an effective potential [fl8| 



V eS (r) = V ext (r) + 2gn 1 (r) + 2g\*(r)f 



(2.2) 



Thermal averages are calculated with a standard Bose 
distribution of the noncondensed particles in chemical 
equilibrium with the condensate governed by the same 
chemical potential p. In particular the density n\{r) is 
given by fflJTJ] 



ni(r) 



1 



(2wh) v J exp[{p 2 / 2m + V eS (r 



p}/k B T]-V 
(2.3) 



where k B is the Boltzmann constan t, a nd T >(= 2, 3) is 
the dimension of space. Equations (2.1) — (2.3) above 
are the principal equations of the present model. The 
total number of particles N of the system is given by the 
number equation 



N = Nn + 



p(E)dE 



exp[(E-p)/k B T] -1' 



(2.4) 



where N = J \^(r)\ 2 dPr is the total number of parti- 
cles in the condensate. The cri tica l temperature To is 
obtained as the solution of Eq. (2.4) with Nq and p set 
equal to 0. The semiclassical density of states p(E) of 
noncondensed particles is given by 



P{E) 



27rm p / 2 
(2irT<) v 



v eS (r)<E 



i8(E-V cS (r))]V- 2 y 2 d v r. 



The average single- particle energy of the noncondensed 
particles is given by 



{E)nc = 



Ep(E)dE 



exp[(E- n)/k B T\ -1' 



(2.6) 



The kinetic energy of the condensate is assumed to be 
negligible and its interaction energy per particle is given 
by (E) c = {g/2) J y 4 (r)d v r. The quantity of experi- 
mental interest is the average energy (E) = [{E) nc (N — 
Nq)/2 + (E)c]/N, which we calculate in the following. 
The specific heat is defined by C = d{E)/dT g. 

Equations (2.1) — (|2.5|) are to be solved iteratively. 



(2.5) 



The iteration is started at a fixed temperature with a 
trial chemical potential p using ni(r) — 0. Then \&(r ) 
and V e a(r) are calculated using Eqs. ( |2.l| ) and (^). 
With these results ni(r), \&(r), and V e a(r) are recalcu- 
lated. This procedure is repeated until desired precision 
is obtained. The results for the lowest order of iteration 
with ni(r) = will be denoted by I = 1, and succes- 
sive orders by I = 2,3,... etc. With the solutions \l/(r) 
and V e Q (r) so obtained, the density of states p{E) of Eq. 
( |2.5D is calculated. Then it is seen if they satisfy the num- 
ber equation ([2.4|). If Eq. (2^) is satisfied the desired 
solution is obtained. If not, the initial trial p is varied 
until the number equation is satisfied. In each order of 
iteration we calculate the condensate fraction Nq/N and 
energy (E). 



III. NUMERICAL STUDY 

First we consider the three-dimensional case. In this 
case the coupling g is given by g = 4:nh 2 a/m, where a is 
the scattering length. Usually, a dimensionless coupling 
is introduced via rj = (mg / 'nh 2 ) / 'a-^ — Aa/a^Q, where 
aj 10 = y^h/muj. The semiclassical model under consid- 
eration should break down as either rj or N is increased. 
This will violate the condition of weak interaction and 
diluteness. In the case of experiment on BEC of 87 Rb, 
r] = 0.025 and N = 40000 To test our calculational 

scheme, first we solve the present model for ?y = 0.025 
and N = 40000. Our results are very similar to those of 
Ref. [|| . The small difference between these two calcula- 
tions is due to the use of an isotropic harmonic oscillator 
potential in this work and an anisotropic potential in Ref. 
P). In this case specific heat is positive at all tempera- 
tures. Next we increase rj and N. We find that as 77 and 
N are increased, the energy develops a maximum at a 
temperature below T . Consequently, the specific heat 
becomes negative for an interval of temperature above 
this maximum. To show this violation in a pronounced 
way here we show the results for the following three cases: 
(a) T] = 0.1, N = 10 6 , (b) T] = 0.1, N = 10 7 , and (c) 
r] = 0.5, N = 10 6 . In the case of 23 Na, the experimental 
N was as high as 10 7 ||. As the interaction is repul- 
sive in both 87 Rb and 23 Na, it should be possible to have 
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10 7 atoms in a BEC of 87 Rb under favorable experimen- 
tal conditions. Hence these values of the parameters are 
within the present experimental scenario. 
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Fig. 1. Condensate fraction Nq/N and energy 
(E) /NksTo in three dimensions as a function of T/Tq 
for (a) T) = 0.1 and iV = 10 6 , (b) r\ = 0.1 and N = 10 7 , 
and (c) -q — 0.5 and N — 10 6 for iterations 7=1 (full 
line), 2 (dashed-dotted line), and 4 (dotted line). The 
straight line represents the classical Maxwell-Boltzman 
result of energy. 

In Figs. 1 (a), (b), and (c) we plot (E)/Nk B T and 
Nq/N versus T/Tq for different iterations for the above 
three cases. We find that the energies are acceptable un- 
der conditions of diluteness and very weak interactions, 
but as 77 and N increases the average energy of the sys- 



tem may have a maximum leading to a negative specific 
heat above the maximum. From Figs. 1(a), (b), and (c) 
we find that this violation happens in the case (c) which 
has the largest 77 and N. We verified that for sufficiently 
large values of rj and the lowest-order energy also leads 
to negative specific heat. 

Next we consider the two-dimensional case. In this 
case, in analogy with the three-dimensional di- 
mensionless coupling is introduced by 77 = (mg /irh 2 ). 
This coupling is already dimensionless, whereas we 
needed to divide it by a nQ in three dimensions to make it 
dimensionless. First we repeat the calculations reported 
in Ref. Q and our results are in agreement with that 
study. In addition, in agreement with our finding in 
three dimensions, we find that for small values of 77 and 
N the equations of the model converge well and lead to 
acceptable values for condensate fraction Nq/N and en- 
ergy (E)/NkBTo. For larger values of 77 and/or N, the 
condensate fraction No/N is quite acceptable with a tem- 
perature dependence similar to that in three dimensions. 
However, the energy produces a maximum as N and 77 
increase. Hence we shall be limited to a consideration of 
energy only. 

As there is no experimental guideline for probable val- 
ues of N and 77 in two dimensions, as in Ref. Q, we 
consider N = 10 5 and 77 = 0.1, 1, and 10. In Fig. 2 
we plot the temperature dependence of average energy 
(E) /NUbTq for different iterations. The lowest-order re- 
sults for 77 = 0.1 and 1 are in agreement with those of 
Ref. The classical Maxwell-Boltzmann result for a 
noninteracting gas is also shown in Fig. 2. For a weakly 
interacting gas with 77 = 0.1, the energies for all orders of 
iteration are acceptable. For a stronger interaction with 
7; = 1, the lowest-order energy is acceptable. However, in 
this case the energies for all orders of iteration produce 
maxima leading to negative specific heat. For 77 = 10 the 
lowest-order result already leads to a negative specific 
heat. 

Fig 2 




Fig. 2. Energy (E) /NksT^ in two dimensions as a 
function of T/T for N = 10 5 and 77 = 0.1, 1, and 10 for 
iterations I = 1,2 and 4. Notations are the same as in 
Fig. 1. The curves are labeled by the values of 77. 

In addition to the trouble discussed above, the mean- 
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field model may exhibit another unacceptable behavior. 
We see in Figs. 1 and 2, that the energy the system 
could be larger than the corresponding classical Maxwell- 
Boltzmann result for a nonintcracting ideal gas. We re- 
call that for BEC to materialize, the energy of the con- 
densate should be smaller than the corresponding energy 
of the noncondensed system. However, one should con- 
sider the correction to the classical result above due to 
the presence of the trap and the interatomic interaction. 
In all our calculations we have noted that UbTq » Tito, 
so that for temperatures close to T considered above 
the correction to the classical energy due to the pres- 
ence of the trap can be neglected. The same is also true 
for the interatomic interaction at temperatures close to 
To, for the values of r\ considered in this study. Conse- 
quently, the energy of the condensate of the mean-field 
model could be larger than the classical result signaling 
a breakdown of the model. 

It is well known that the semiclassical model requires 
ksTo >> hu> and the Thomas- Fermi approximation does 
not hold for very small number of atoms in the con- 
densate. We have verified that for the cases studied 
here /c^To is typically 50 to 100 times larger than hu> 
and Thomas-Fermi approximation apparently seems to 
be valid for the condensates of size N ~ 10 6 to 10 7 con- 
sidered here. The diluteness condition for the validity of 
the GP equation na 3 << 1 is also valid in three dimen- 
sions, where n is the density. For the case considered in 
Fig. 1(a), N = 10 6 , a = 50 nm, for a condensate of typ- 
ical dimension 100 jtim, na 3 ~ 10~ 3 << 1. The values 
of na 3 are larger for Figs. 1(b) and 1(c). Nevertheless, 
we find that these values of diluteness and coupling set a 
limit to the applicability of the semiclassical mean-field 
two-fluid models for studying BEC thermodynamics. 

IV. CONCLUSION 

In conclusion, we reexamined the problem of BEC un- 
der the action of a trap potential using a two-fluid mean- 
field model H[| in both two and three dimensions. We 
employed an iterative solution scheme of the system of 
equations. Although the system leads to rapid conver- 
gence for a weakly interacting dilute system, with the 
increase of coupling and particle number the iterative 
scheme leads to physically unacceptable results for ther- 
modynamic observables. Specifically, this may lead to 
a maximum in energy responsible for a negative specific 
heat. In three dimensions the breakdown of the mean- 
field model happens for values of coupling and particle 
number, which are not so remote from present experi- 
mental scenario for repulsive interatomic interaction. In 
addition, the energy of the condensate could be larger 
than the corresponding classical energy of the system, 
which is another independent unacceptable result of the 
model. The larger values of the coupling r/ and N consid- 
ered in this work possibly sets a limit to the applicability 



of the mean-field equations for BEC. Summarizing, the 
most important finding of this study is that the itera- 
tive solution of the mean-field model of Refs. ||,[| for 
a condensate may lead to unphysical thermodynamical 
properties for medium to large coupling and number of 
particles. 

Despite the above deficiency of the mean-field two-fluid 
thermodynamical models, they continue to be very use- 
ful in many cases. The virtue of these models is the 
simplicity and ability to yield results in agreement with 
experiment for weak interactions and dilute systems. As 
a theoretically sound description of the BEC thermo- 
dynamics seems to be unmanagably complicated, these 
simple mean-field two-fluid models remain as attractive 
simple alternatives to study thermodynamical properties 
provided that proper attention is paid to remain inside 
the domain of their validity. 
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